Modeling Bose Einstein Correlations via Elementary Emitting Cells 
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We propose a method of numerical modeling Bose Einstein Correlations by using the notion of 
Elementary Emitting Cells (EEC). They are intermediary objects containing identical bosons and 
are supposed to be produced independently during the hadronization process. Only bosons in EEC, 
which represents a single quantum state here, are subjected to the effects of Bose-Einstein (BE) 
statistics, which forces them to follow a geometrical distribution. There are no such effects between 
particles from different EECs. We illustrate our proposition by calculating a representative number 
of typical distributions and discussing their sensitivity to EECs and their characteristics. 
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I. INTRODUCTION 

In every high energy colhsion experiment, a vast num- 
ber of secondaries (mostly pions) is produced encoding 
valuable information on the dynamics of the hadroniza- 
tion process. Because of their complexity, such reactions 
can only be investigated by numerical modeling using 
specialized codes, Monte Carlo event generators (MCEG) 
They all use positively defined probabilities to de- 
scribe the particle production process. This means that 
they neglect ofF-diagonal terms in the corresponding den- 
sity matrices formally describing such processes. As a 
result they did not properly implement all the correla- 
tions, in particular the quantum mechanical correlations 
between identical particles (both of Bose-Einstein type 
for identical bosons (BE) and of Fermi-Dirac type for 
identical fermions (FD)). The whole collision process is 
expressed only by the product of single particle distribu- 
tions. This makes the present MCEG a priori impossi- 
ble to properly account for results of many experiments 
(starting from [3]), which explicitly show features usually 
attributed to correlations of this type 0]. To describe 
such experiments, one has to introduce, in some way, ef- 
fects of quantum statistics, opening thus a new subject - 
Bose-Einstein (BEG) or Fermi-Dirac correlations Q. 

This subject already has a long and well documented 
history 0] but it is still far from complete. Its most 
specific features one is faced with are the following: In 
experiments of an exclusive type (measuring all produced 
secondaries, as in 0]) it is possible (at least in principle) 
to construct the properly symmetrized wave function of 
all identical pions and thus account for our inability to 
determine which particle is emitted from which position 
in the source. In the simplest approach using a plane 
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waves representation of particleSj a such wave function 
for riir-pion state, has the form 
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E 



exp 



(1) 



{<j{j) denotes the j*'' element of permutation of the se- 
quence {1, 2, n.^.} and we sum over all n^! permuta- 
tions in this sequence, including the identity permuta- 
tion). Points of detection, x's, will vanish when calculat- 
ing probabilities, but points of production of secondaries, 
r's, will remain. Because experiments measure only mo- 
menta of particles, one must integrate over {^j} using 
some multiparticle spatio-temporal distribution function, 
p{{rj}). This is usually assumed to be factorizable, i.e., 
expressed by the product of independent single particle 
distributions, p{{rj}) — Yij Pi''^j) H- Assuming further 
a totally chaotic hadronizing source one gets the proba- 
bility of the n^-pion state in the form of permanent of 
the matrix 



'Pi, 
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with 



''->(r)dHr}, q,,^p,-p„ (3) 



(which are all real if p(r) = p(— r)). 

However, most of modern experiments are of an in- 
clusive type as one measures effect of BEG on limited 
samples of secondaries only. The unobserved part of 
the system acts then as a kind of thermal bath influenc- 
ing measured samples of data [§]. Also the registered 
multiplicities are much higher prohibiting calculation of 
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permanents given by eq. H]). The subsystem subjected 
to further experimental (and theoretical) analysis forms 
only a (not too well defined) part of the total system 
with strongly fluctuating characteristics to be averaged 
during data collection procedure. Because, as we have 
mentioned, such systems can be described only by nu- 
merical modeling, one faces the fundamental question of 
how to account for BEC in such modeling processes? It 
must be realized that this is a completely different thing 
then a simple calculation of the usually considered n^r = 2 
case, which in the case considered here amounts to a well 
known expression for the 2— particle BE correlation func- 
tion defined as a ratio of a two particle distribution and 
single particle distributions with Q = \pi — Pj\ , 



C2{Q) = 
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C^—^ 2 when Q and C2 ^ 1 for large values of Q) 
0, B @|- In fact, in the majority of cases, one is fol- 
lowing precisely this route with different functional forms 
of p{QR) and sometimes even with an entirely different 
form of C2 (cf., for example, fio'l), but always derived in 
some analytical way - both for one and three dimensional 
cases. One is then fitting such correlation functions to 
appropriate sets of data, aiming to get information on 
the quantity R (and similar others). The reason for this 
is that C2{Q) is usually regarded as a (kind of) Fourier 
transform of the space-time characteristic of the emitting 
source, p, therefore providing information on p(r) jllj . 

It is precisely this fact that makes BEC so interesting, 
but we shall not follow this route here. Our aim is to 
provide an algorithm that could cope with BEC from its 
first steps (to model it) and which could bring this ef- 
fect to other MCEGs once implemented there. This we 
shall do by using the notion of Elementary Emitting Cells 
(EEC) introduced some time ago [T5| . They are some in- 
termediary objects (in fact quantum states) containing 
identical bosons (henceforth we shall assume them being 
pions), which are supposed to be produced independently 
during the hadronization process. Only bosons in EEC 
are subjected to Bose-Einstein (BE) statistics, which 
makes them follow a geometrical distribution. There are 
no such effects between particles from different EECs. 
Using this idea we propose an algorithm that can be 
used in any MCEG, which has effects of BEC as one 
of its basic features. This will be presented on the back- 
ground of previous attempts to model BEC discussed in 
[I, E Hi, In, [III- We shall illustrate its action by cal- 
culating a representative number of typical distributions 
and discuss their sensitivity to EECs and their charac- 
teristics. In this work no comparison with experimental 
data is offered, we limit ourselves only to a thorough dis- 
cussion of our algorithm. It is supposed to be the main 
building block of any serious MCEG aiming for a real 
comparison with data and including therefore, for ex- 
ample, also the distribution of energies one is supposed 



to hadronize or possible flows in the system, i.e., sub- 
jects which are out of the scope of this work. In particu- 
lar it will be most useful in cases where the hadronizing 
source is well deflned (as in e+e^ annihilation processes 
and other elementary processes) or where the number 
of secondaries is exceptionally large (in high multiplicity 
events) . 

For completeness, in the next Section (|TT| we shall pro- 
vide a short overview of numerical modelling of BEC pro- 
posed so far. This will put our investigation in proper 
perspective. Section Hill contains details of the proposed 
algorithm and examples of results obtained from its sim- 
plest version. The last Section contains conclusions and 
a discussion. Some speciflc problems connecting with the 
method proposed are addressed in Appendices [Xl- [Cl 



II. HISTORY OF NUMERICAL MODELING OF 

BEC 

A. Imitating BEC with existing MCEGs 

Imitating BEC means to use the original outputs of 
the existing MCEG codes and change them in such way 
as to reproduce measured C2{Q) (or other experimen- 
tal characteristics of BEC). Two approaches are used for 
this. In the first one one introduces some special global 
weights (i.e., weights built for the whole event) to bias 
accordingly the original results of MCEG [l3| (usually 
checking whether other observables were not changed too 
much - otherwise one has to re-run MCEG with new in- 
put parameters). This procedure is justified by noticing 
the following [1^ : Let M — M^, be the matrix ele- 
ment describing the production of a hadronic final state 
of n identical bosons. It consists of n\ terms, each corre- 
sponding to a particular permutation a of the n identical 
particles in the final state. In the simulation process of 
MCEG the interference terms (off-diagonal elements in 
permanent eq. ([2])) are neglected and one gets that the 
probability to produce such a state is 



\M\ 



MCEG 



(5) 



To remedy this situation one assumes some weight 
assigned to each event such that 



(6) 



There is no unique way to choose the weights Wa- The 
only requirements is that one gets good fits to the corre- 
sponding C2 functions [13, [ll] ■ The tacit understanding 
is that then also ~ |Mp. 

In the second approach, one locally modifies (by 
weighting each pair of particles in a given event) the orig- 
inal output of the MCEG used. This can be done either 
by modifying its energy momentum spectra Il9| or by 
changing the resulting charge assignment [l^, l2ll |. In the 
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first case one introduces weight function fsEiQ) for pair 
of particles momenta which are changed, such that 



/•Q rQ+SQ 
/ dn{q) = / fBE{q)dniq), 
JO Jq 



(7) 



i.e., for fBsiq) > one has SQ < 0. In this way 
the energy-momentum imbalance that results from such 
a procedure is properly accounted for and number of 
particles, N = J dH., is conserved [l^. In the second 
case [lO, HH the original spatio-temporal and energy- 
momentum structure of the original event is preserved, 
but often spurious unlike particle correlations occur [2^ . 
What must be stressed is that this approach, contrary to 
the previous one, always works on the level of a single 
event. It is therefore more suitable as an additional tool 
(sometimes called afterburner) to be used together with 
the known MCEGs. 

It must be stressed that all these methods modify 
original physics underlying MCEGs in an essentially un- 
known way (23j . Using them one assumes therefore that 
changes incurred are small and irrelevant. We shall pro- 
ceed now to attempts to simulate BEG by which we un- 
derstand a situation in which the algorithm introducing 
effects of quan tum statistics involves all produced parti- 
cles 



ts or quan tum 



B. Attempts to simulate BEC numerically 

1. Metropolis importance sampling mettiod 

In the standard Monte Garlo technique due to 
Metropolis [Metropolis importance sampling method) was 
used. This is a general method to generate an ensemble 
of n-body configurations according to some prescribed 
probability density. In 8] this technique was used to 
modify the directions of momentum vectors of selected 
particles from a system of n identical particles in or- 
der to impose the n-particle distributions derived from 
BE correlation functions. In particular, it was done 
by changing the momenta of selected particles, pi — > 
p'j €: d'^N/dp^, in such way as to maximize the probabil- 
ity of detection of the n^-multiparticle state, 'Pi,...,p„ , 
i.e. accepting a new configuration with probability P = 
min{l,Pnew/Poid}, where Pom = V{1, ...,Pi, and 
= V{1, ...,p'i, ...,Pn„}- This procedure is then 
repeated many times until a kind of "equilibrium" is 
achieved. As shown in [sj , one was able in this way to gen- 
erate typical multipion events, which explicitly exhibit all 
correlations induced by Bose statistics. The most impor- 
tant result for our further consideration is the fact that 
as a result of the application of this algorithm a number 
of objects, called speckles in Q and being clusters of a 
number of identical pions in phase-space, is formed. It 
means that in the multidimensional phase space perma- 
nent ([2]) exhibits rich structure of local maxima (attract- 
ing particles) and voids (repulsing them), which replaces 



the original distribution one started from. Actually the 
only drawback of this method is that symmetrization of 
clusters with sizes larger than Uduster ~ 10 takes a pro- 
hibitively long time. 

Two points must be stressed when summarizing this 
symmetrization procedure. The first is that it involves 
all (identical) secondaries in the event under considera- 
tion, some producing specific structure in their distribu- 
tion in the allowed phase space, namely it is clustering 
them in some regions of phase space. The second point is 
that this phenomenon leads immediately to a broadening 
of the resultant multiplicity distribution (MD): starting 
from a Poisson MD for a non-symmetrized wave function 
one ends up after symmetrization with a geometrical (or 
Bose-Einstein) MD for a single speckle and with Negative 
Binomial MD [25] for the whole system. 

We close this section with the following remark. So far 
particles were represented, for simplicity, by plane waves. 
However, this approach leads to some unpleasant effects 
because it violates the Heisenberg uncertainty relation 
constraining the simultaneous specification of coordi- 
nates and momentum as implied by eq.([4]). For example, 
in the case of sources with strong position-momentum 
correlations the two-particle correlation function C2{Q) 
can drop significantly below unit y [1 , [26j . This method 
has therefore been generalized in [13|, where plane waves 
have been replaced by wave packets. Both features men- 
tioned above were also observed here. Therefore in 
modification simplifying the selection process was pro- 
posed. It was argued that it could be limited only to 
identical particles whose wave functions overlap in phase 
space, i.e., to particles forming speckles or clusters men- 
tioned above (with the size of this overlap being a new pa- 
rameter) . Notice that such decomposition corresponds to 
replacing the full permanent in eq. ^ by matrix with a 
block structure, each block representing one cluster with 
no correlations between them: 



EECi 











EEC 



(8) 



In what follows we shall identify these blocks with EECs 
mentioned before. However, it should be remembered 
that cells selected here were so far preselected in ($,9) 
space only and, by construction, tend to contain parti- 
cles with similar momenta. This method differs therefore 
from that we are going to describe now in which EECs 
are created dynamically without any restriction p^ . 



2. The acceptance-rejectton method 

The acceptance-rejection method used in (isj is the 
well known " hit-or-miss" technique of generating a set of 
random numbers according to a prescribed distribution, 
here given by an expression describing the collapse of a 
multiparticle wave function into a properly symmetrized 
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state represented by eq.(l2]), as required by Bose quan- 
tum statistics. In contrast to that discussed above this 
method is sequential because the n multiparticle event 
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FIG. 1: Schematic example of acceptance-rejection algo- 
rithm proposed in [l^ at work; Ci+i(max) — (i -I- 1)! and 
R £ [0, 1] is random number. For i particles already present 
with Pi^i fixed (gray circles, here i = 9) one calculates 
Ci+i({pi,...,i},Pa: = Pi+i) and selects a random number R. 
For the situation shown here (with schematic form of C(px) 
as example) a new particle can only be added to "cells" 1, 3 
and 5 but not to 2 and 4. 

is constructed by adding the k*^ particle to the (fc — 1) 
particles already selected, until fc = n is reached. This is 
done in the following way. One starts with empty phase 
space and inserts in it the first particle with momentum 
chosen according to some distribution (for example, the 
one which is supposed to reproduce single particle spec- 
tra). The presence of the first particle causes that a sec- 
ond particle must be added according to the 2-particle 
distribution function, which for identical bosons is given 
by the 2-particle correlation function C2 oc P12 with mo- 
mentum of the first particle fixed and momentum of the 
second particle to be selected. Because C2{Q — Pi — P2) 
has a maximum at Q = 0, the second particle will tend 
to be located near the first one in phase space but a pri- 
ori it can take any position in it. Addition of a third 
particle must be now performed according to C3 (x Pi, 2, 3 
(with momenta pi and p2 already fixed and only ps being 
selected), which again has bumpy structure, especially 
when particles 1 and 2 are located far from each other. 
Therefore the third particle will most probably locate it- 
self near one of the already selected particles but a priori 
it can take any position in phase space. If, by chance, 
it will be far away from both particles 1 and 2, it will 
become for future particles a new point of attraction and 
seed of the new bump. Notice that if, say k particles oc- 
cupy the same region of phase space, the strength of the 
bump they form, which attracts other particles, increases 



fc! times. This process, schematically illustrated in Fig. 
[TJ continues until all particles are used (their number can 
be either preselected, in which case the initial energy will 
vary, or can result from the procedure itself when initial 
energy is fixed). 

To summarize: this procedure leads again to some spe- 
cific nonuniform distribution of particles in the allowed 
phase space with some cell-like structures (bumps) show- 
ing up. It results from the fact that regions with some 
particles inside them already present will have a bigger 
chance to attract a new particle. Unfortunately, this se- 
quential procedure is even more time consuming than the 
previous one and therefore rather unpractical. 



3. Information Theory approach 

The only workable example of MCEG with features of 
Bose-Einstein statistics built in from the very beginning 
has been proposed in [16]. The known total energy Etot 
has there been distributed among a given (mean) num- 
ber of secondaries, n ~ f;,_|_ + n_ + uq (where = "■-), 
with limited transverse momentum parameterized by the 
mean value (pt), and it was done in such way as to re- 
produce data on both single particle distributions and 
those for BEG as well. To this end Information The- 
ory (IT) method was used to obtain the most proba- 
ble (and least biased) formula for rapidity distributions, 
dN/dy, and multiplicity distributions P{n). It resembles 
the usual grand canonical distribution but is more gen- 
eral than that because the "temperature" T and "chem- 
ical potential" /i occurring there are now two lagrange 
multipliers obtained by solving the corresponding en- 
ergy and particle conservation constraints. To also get 
C2{Q) it turned out that it is enough to additionally di- 
vide the available rapidity space into cells of fixed size 
5y each and assumption that each cell can contain only 
identical particles (i.e., pions of the same charge in this 
case). It is remarkable that, in addition to reproducing 
all single-particle characteristics of the collision as well as 
backward-forward correlations, with this step one gets at 
the same time also multiplicity distributions P{n) in Neg- 
ative Binomial (NB) form [23], the proper structure of 
the twojgarticle BEG function C2{Q) and intermittency 
signal 28] . The distinctive feature of this method is that 
it deals only with measured momenta of produced secon- 
daries, there is no trace on unmeasured spatio-temporal 
structure present in other methods. If one now wants to 
somehow deduce this information one has to treat C2{Q) 
in the same way one is treating all experimental results 
on BEG (i.e., in fact one has to assume that it is a kind 
of Fourier transform of the hadronizin g so urce p{r) and 
perform its routine analysis as in 0, 

The most important result of [l6| for us is demonstra- 
tion that the decisive factor leading to proper structure 
of the correlation function C2{Q) was bunching property 
in the rapidity space introduced there. 



5 



C. Quantum optical analogy 

Let us therefore concentrate for a moment on this 
bunching feature of bosons introduced in [l^. At first 
let us remind ourselves that it has been noticed and dis- 
cussed already many times HO, HU, IHl and was re- 



garded as manifestation of quantum statistics 3J| . It is 
especially widely discussed and used in quantum optics 
[33] but it was also employed to describe some aspects of 
the multiparticle production processes [H, [H, Ha] • This 
is because C2 can be also regarded as some measure of 
correlation of fluctuations. One uses here the fact that 

(nin2) = (ni)(ri2) + (("-1 - ("-1)) ("-2 - ('^2))) 

= {ni){n2) + pa{ni)a{n2) (9) 

(where a{n) is dispersion of the multiplicity distribu- 
tion P{n) and p is the correlation coefficient depend- 
ing on the type of particles produced: p — -1-1,-1,0 for 
bosons, fermions and Boltzmann statistics, respectively). 
It means therefore that one can write the two-particle 
correlation function ([4]) in terms of the above covariances 
© stressing therefore its stochastic character (20l. [2TI. [37| : 



C2{Q) 



(n,; iPi))(nj (pj)) 
a{ni) cr{nj) 



P 



[rii ipi)) {Uj (pj)) {n^{p^))' 



(10) 



Notice that for geometrical (Bose-Einstein) multiplicity 
distribution (rii) = {nt (pi) [rii (pi) — 1]) (correspond- 
ing to bosons, p = 1) one gets, for i — j, C2{Q = 0) = 2, 
i.e., its maximal allowed value. This bunching prop- 
erty has already been used in our previous proposition of 
MCEGof the "afterburner" type [13, [HI mentioned in 
Sec. | II Al it was realized in the form of EECs introduced 
in |l2|. Essentially the same idea will be the cornerstone 
of our algorithm, which we shall now present. Identi- 
cal pions will be assumed to be subjected to BEG only 
when inside EEC, those from different EECs are totally 
independent (this can be also expressed using notation 
chaotic and coherent, see Appendix |^ . 



III. MODELING BOSE EINSTEIN 
CORRELATIONS VIA ELEMENTARY 
EMITTING CELLS - OUR ALGORITHM 

A. Introduction 

The lesson learned from the approaches presented in 
Sec. |TT] is that construction of the proper quantum mul- 
tiparticle bosonic state can be performed (i) either by 
symmetrization of the corresponding wave function con- 
structed for all produced particles @, [H, [13, [3| or (ii) 
by following quantum statistical arguments formulated in 
[sol . I31I [32 :] and directly bunching produced identical sec- 
ondaries with (almost) the same energy in phase space to 



form EECs [H, [13, [il| . So far, the fact that the distribu- 
tion obtained this way is of NBD type was only shown in 
[1] . Nevertheless the emergence of bunching (both in ^ 
and [isj ) is convincing and makes it an essential charac- 
teristic of the bosonic character of produced secondaries 
one has to account for. This will be our main assumption 
in what follows. 

We have decided to model the effects of proper sym- 
metrization of a multiparticle state by assuring that iden- 
tical particles are produced in bunches with geometrical 
distribution of particles. To get such a distribution one 
has to choose particles sequentially with some prescribed 
probability V until first failure. If this failure happens 
for the {N + 1)*'^ trial one gets immediately that 



P{N) = (1 - V)V'^ with (N) 



l-V 



Notice that for V defined as 



V = V[j ■ exp 



El 

' T 



(11) 



(12) 



(and only then) one gets the usual form of the Bose- 
Einstein distribution for the i*'' EEC: 



{N{E,)) 



1 



(13) 



Because, according to eq. (|T2t Vq controls the maximal 
number of particles which can be allocated in a given 
EEC, one can formally introduce a kind of "chemical po- 
tential" (as in [1^) defined as iJL — T\nVo and get j38|] 



1 



1 



(14) 



Such a choice of V is therefore very crucial in the process 
of further modeling BEC and is thus the cornerstone of 
our algorithm. 



C,(n) 



(a) 



C,(Q) 



(b) 
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C,(Q) 






(c) 













5 10 15 20 25 30 35 40 45 50 1 2 3 4 5 5 10 IS 20 25 30 

n ° [GeV] Q [GeV] 

FIG. 2: Example of C2 occurring for a picnic lattice in (one- 
dimensional) momentum space: (a) C2(n) as function of num- 
ber of particles in a given cell for different A'^; (b) C2{Q) as 
function of Q (in both cases V — 0.5). In (c) C2(Q) as func- 
tion of Q defined by eq. (|15|l for artificially limited cell occu- 
pancy, i < Umax (now V — 0.9). 

Let us now demonstrate this procedure on a simple 
example of a one-dimensional pionic lattice model [2l|. 
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The pions located in the sites of the lattice (which pre- 
vents identical pions to be found at the same point of 
phase-space) are endowed with charges in such way that 
after selecting the charge of the first pion, one assigns the 
same charge to the consecutive neighboring pions with 
probability V until first failure. After failure one assigns 
(again randomly chosen) charge to the next pion in line 
and continues this procedure until all pions are used. In 
this way a number of cells is formed with particles in 
them distributed by construction (cf. eq. [TT|) and ^,38]) in 
geometrical fashion, cf., Fig. [51 In Fig. dependence 
of C2 on the lattice spacing defined by different values of 
N is presented whereas Fig. [2Jd presents its dependence 
on Q defined as 

Q=\Pi- Pj \= ■\i-j\=Sp-n (15) 

(where p e [-Pmax^Pmax] with pmax = 10 GeV and 
= 100 denotes the total number of sites in our lat- 
tice) . Comparing them one can deduce that the width of 
the correlation function, (t(C2), is roughly proportional 
to the product of the average number of particles in the 
cell, (Npart), and the average distance (Sp) between these 
particle on the lattice, cr(C2) c>c (Npart) ■ {^p) ■ Because the 
average distance between particles on the lattice is fixed 
by the (constant) lattice spacing, the width of the corre- 
lation function depends only on the mean cell occupancy, 
(Npart), which is directly related to the probability V (cf. 
eq. (fT3|l ) and for constant V all widths of C2{n) in Fig. [2k 
are the same. Finally, in Fig. [2}: the correlation function 
C2{Q) is shown for different limits Nmax imposed on the 
maximal allowed cell occupancy (i < Nmax) (for V = 0.9 
to assure large occupancy of EECs). Notice that both 
the value of the intercept parameter, A = C2((5 = 0) — 1, 
and the width of the correlation function C2 depend on 
the maximally allowed number of pions in one cell, Nmax- 



B. Implementation 

This experience prompted us to propose that BEC 
should be introduced into the modeling procedure of mul- 
tiparticle production processes as early as possible, ide- 
ally at the very first stage of the selection procedure. It 
means that one must devise some procedure how to di- 
vide a given amount of energy W into produced bosonic 
particles, assumed to be pions of charges: (-)-, — , 0), with- 
out any a priori assumption on what concerns their mul- 
tiplicities other then n*^+) = but with effects of quan- 
tum statistics accounted for. Among procedures men- 
tioned above only that in Sect. Ill B 21 [Ti| satisfies this 
demand, however, because it always involves all produced 
particles, it very quickly becomes impossible to follow be- 
cause of CPU time demand. On the other hand, results 
of [1, [1^ show explicitly that quantum statistics leads to 
bunching particles in phase space. Such bunching was 
therefore assumed in [l^, [ij, [la| and in fT(l| , for the first 



time, the geometrical distribution form of particles in 
bunches was used p7| . 




FIG. 3: Schematic view of Quantum Clan Model (QCM). 

We shall assume therefore that particles are produced 
according to a mechanism which can be regarded as a 
quantum version of the clan model (CM) introduced in 
[25| to explain the fact that apparently all multiparti- 
cle production processes result in the NB form of result- 
ing multiplicity distributions P{N); we call it therefore 
the Quantum Clan Model (QCM), cf. Fig. [31 In CM 
model (distinguishable) particles are supposed to be pro- 
duced in clans, which were supposed to be formed inde- 
pendently (and therefore distributed in Poisson fashion), 
with particles in clans following a logarithmic distribu- 
tion. Convolution of these two results in the NB form of 
P{N) {V denotes the probability of producing a particle), 

PNB{N)^[^^l^^^Y\l-Vr. (16) 

In QCM, because of quantum statistics, we must assume 
that each quantum clan represents a single quantum 
state and therefore contains only identical particles of 
(almost) the same energy, which are distributed geomet- 
rically [H, [1^ [S^l • Keeping the assumption of indepen- 
dent production of quantum clans, one immediately finds 
that previous NB multiplicity distribution (|16p is now re- 
placed by the so called Polya-Aeppli (Geometric-Poisson) 
distribution defined as §^ (with Q ^ {I -V)(N)): 

Ppa{N) = Poisson Geometrical (17) 

which was already used in multiparticle phenomenology 
some time ago w]. It differs from the NB only at small 
multiplicities. Otherwise it is essentially the same. 

To implement numerically the QCM we proceed in the 
following way. 

(i) A pion of randomly chosen charge is selected with 
energy Ei — Eeec following some assumed distribu- 
tion f(E). The form of this distribution is dictated by 
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FIG. 4: Examples of distribution of particles in a given EEC 
(upper left panel), distribution of EECs (upper right panel) 
and total charged particle distribution resulted from the con- 
volution of distributions shown in the first two panels (lower 
panel, open circles and solid curve, see text for details). No- 
tice that, with this choice of parameters, one can reproduce 
exactly the experimental data for e~e"'" annihilations taken at 
the same energy [JJ] (black circles). 



physics of the model used to describe the production pro- 
cess which we would like to follow (much the same as 
in [IBl)- It has no influence on the bosonization (other 
then emerging from conservation laws), on the other hand 
it strongly affects the final form of multiparticle distri- 
butions obtained (by changing effective distributions of 
EECs). This pion is therefore supposed to be the seed 
for the first EEC. 

(a) The other pions of the same charge are than added, 
one after another, using probability V as given by eq. 
P^ . This is done until the first failure, which marks the 
end of formation of the first EEC. After that, one starts 
formation of another EEC by proceeding to point (i) and 
choosing from f{E) another energy Eeec of another pion 
of randomly chosen charge, which forms the seed of the 
second EEC. The whole procedure is repeated until all 
available energy M is used up. 

{Hi) Once accepted, each of the selected pions forming 
the first EEC is endowed with energy Ei {i > 2), which 
is selected from some distribution G{Ei) centered on the 
energy of the first pion forming the seed of this EEC, 
Eeec- The sense of the function G{Ei) is that it reflects 
the fact that the requirement that all particles belong- 



ing to a given EEC must be in the same energy state, 
in which case one expects that G{Ei) ~ 5 [Eeec — Ei), 
means that their energies can differ only by an amount 
corresponding to the width of the spectral line mentioned 
in [31] . We shall then use G{E) either in the form of delta 
function mentioned above or parameterize it by a Gaus- 
sian, 



G{Ei) = 



1 



exp 



{Eeec — Ej) 
2a2 



(18) 



In the thermal- like model for f{E) used here, there is a 
natural length scale given by the temperature T and it is 
therefore natural to choose <t as being proportional to it, 
we shall therefore take a = <jqT (for the possible physical 
meaning of a see Appendix [B]) . 



TABLE I: The mean total multiplicities and their dispersions 
((A^) and (Jn), the mean total multiplicities of EECs and their 
dispersions {{Ncm) and ctn^^u) and mean multiplicities and 
dispersions in EEC {{Npart) and UNpart ) calculated for some 
selected choices of parameters Vo and T (upper part) and for 
different initial energies {W{i) = 0.5W, W{ii) = W = 91.2 
GeV and W{iii) = 2W - lower part). All results are for ao = 
0.0 except the middle part, which was obtained for ao — 0.3 
for comparison. 
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0.9 




41.39 


12.87 


18.21 
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2.27 


2.61 


0.7 


3.5 


33.81 


8.28 


20.58 
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1.28 


0.5 




30.34 


6.63 


22.41 


4.11 


1.35 


0.78 


0.3 




28.22 


5.69 


24.01 


4.44 


1.17 


0.48 
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82.68 


15.29 


39.85 


4.66 


2.07 


2.13 
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3.5 


41.39 


12.87 


18.21 


3.24 


2.27 


2.61 




5.5 


27.99 


11.13 


12.03 


2.64 


2.32 


2.78 


0.7 


3.5 


32.05 


7.12 


19.54 


3.47 


1.63 


1.27 
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(TV) 
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(A^cen) 




{Npart) 






To 














{i) 


3.5 


21.22 


9.08 


9.49 


2.29 


2.24 


2.56 


{ii) 


0.9 


41.39 


12.87 


18.21 


3.24 


2.27 


2.61 


{Hi) 




81.70 


18.27 


35.62 


4.58 


2.29 


2.65 



In this work we shall use the function f{E) in the form 
of a Boltzmann distribution, 



/(i?)=exp(-| 



(19) 



corresponding to a kind of thermal-like model with tem- 
perature T being the main parameter. In this case the pa- 
rameter Vq, governing the number of particles in EECs, 
plays role of the chemical potential. As expected, one 
gets then EECs distributed according to a Poisson distri- 
bution, each EEC containing particles of the same charge 
only, which are distributed according to Bose-Einstein (or 
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geometrical) distribution, also the shape of charged par- 
ticle multiplicity distributions comes out as expected, see 
Fig. m It was obtained for hadronizing energy W = 91.2 
GeV using T = 3.5 GeV, Vo = 0.7 and cto = 0.3. For 
the corresponding values of mean multiplicities and their 
dispersion cf. Table [J (the middle row). 

The W = 91.2 GeV value of energy will be used 
throughout the paper (chosen to allow for the only com- 
parison with data we show in Fig. U}. In all examples 
shown here the number of MC trials was Nmc = 50000 
and reference frame used to calculate C2 function is com- 
posed from (H — ) particles whereas C2 themselves are 
calculated for ( ) pairs. Table |T] shows the correspond- 
ing multiplicities (and their dispersions) of all particles 
obtained when hadronizing mass W and those for the to- 
tal number of EECs and particles in them. The bigger 
Vo, the bigger the total multiplicity, smaller the number 
of EECs and bigger their occupancy. Increase of T de- 
creases the total multiplicity and number of EECs but 
slightly increases their occupancy. The increase of avail- 
able energy W results in an increase of all these quanti- 
ties, except for the cell occupancy, which remains essen- 
tially the same. 




E. [GeV] E, [GeV] 



FIG. 5: Comparison of energy distributions obtained using 
zero and nonzero values of ctq compared (dashed-line) with 
the corresponding Bose-Einstein form of energy dependence 
of occupation number as given by eq. (|13p . The mass of 
hadronizing source is W — 91.2 GeV, except the last curve 
where it is 4-fold greater for comparison. 

Fig. O shows an example of the corresponding energy 
distributions of produced secondaries. Notice that the 
effect of bunching (i.e., effect of introducing EECs) is 
visible only in the limited range of the allowed phase 
space, concentrated at small energies. In the case consid- 
ered here, where the allowed range is (1/2) • {W = 91.2) 
GeV, it practically vanishes for E > 7.0 GeV and af- 
ter that value the distribution follows the exponential 
form of f{E) we have started from. This means that 
BEG increases the multiplicity of event by adding parti- 
cles with small energies (see also Table HI]) . Notice that 
for nonzero uo one gets a displaced maximum for small 
values of E. At large energies results follow the shape of 
the original f{E) distribution used here. The other 



important feature is the fact that none of the numerical 
simulations reproduces the Bose-Einstein form of energy 
dependence of occupation number, usually used in all 
analytical estimations and given by eq. (|13p . This is 
because of the finitcness of available energy W one can 
use for hadronization, which results in limited occupancy 
of EECs and violates conditions used to obtain eq. (fT5|) 
[3^ . Detailed results on the mean number of EECs and 
charged particles multiplicity in them in different energy 
bins are presented in Table [III 

TABLE II: The mean number of EECs, Nceii, and charged 
particles multiplicity in them, n^- , in different energy bins 
calculated for T = 3.5 GeV and two sets of (cro,'Po). 



Bins in E 


(TO = 0.1; Vo = 0.9 


(70 = 0.3; Vo = 0.7 


(GeV) 












Ncell 




Ncell 




0.0-0.5 


0.58 


2.62 


0.63 


1.2 


0.5 - 1.0 


0.72 


3.00 


0.79 


1.60 


1.0-2.0 


1.16 


2.98 


1.26 


2.62 


2.0-3.0 


0.87 


1.58 


0.95 


1.70 


3.0-5.0 


1.15 


1.64 


1.25 


1.76 


5.0 - 7.0 


0.65 


0.77 


0.71 


0.83 


> 7.0 


0.84 


0.89 


0.92 


0.97 




\ [GeV] 6, [GeV] 



FIG. 6: Illustration of importance of spreading in energy. 
Left panel: C2{Se) case with of ao = 0.0, 0.01 and 0.1. For 
(70 = the maximum of C2 is divergent (all points are in 
the first bin). Right panel: C2{Se,p^) cases with ao — 0.0 
or (70 = 0.1. Notice that C2{Sp^) (calculated for momenta 
distributed isotropically) has nonzero width even for ao — 0. 
Introducing nonzero ao results in further broadening of C2. 

We now proceed to the correlation functions C2{Sx — 
Xi — X2) and start with distributions in energy, C2{Se) 
. In Fig. [6] one observes that when using uo = (i.e., for 
strictly (5-like form of G{E)) the whole effect is located in 
the first bin only (this is just computer realization of delta 
function). Therefore, if nonzero widths of C2 are needed, 
one must use ctq > 0. Notice, however, that this width is 
not equal to the input a — a^T used because the differ- 
ence of two variables, each following the same Gaussian 
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distribution, is again Gaussian but with twice cr^, there- 
fore the final distribution should be ~ ^72 broader, which 
is roughly the case shown in Fig. [S] (where a — 0.35 GeV 
was used as input in (jlSp whereas the width which can 
be read off from the obtained shape of C2 is equal to 0.45 
GeV). 

We must stop for a moment to comment the fact that, 
in all figures presenting C^i^x) shown here, their values 
are significantly smaller than unity for large values of Sx ■ 
This is not an artifact of our algorithm, but results from 
the method of presentation of our output. We want to 
keep the same number of pairs both in the real event and 
in the reference one (which, in our work, is always built 
from pairs of opposite charges). Technically it means that 
one has to conserve the area under each curve for C2{Sx)- 
Actually this effect is also known in all other approaches 
to BEC and is usually corrected by arbitrarily shifting C2 
in such way that it equals unity for some large value of 
the argument (in practice set to be equal 1 GeV) [13] ■ We 
shall not do this here, but this fact must be remembered 
when looking at our results. 

In Fig. MC2{Se) is compared with C2{Sp^) calculated 
assuming an isotropic distribution of momenta of par- 
ticles in a given EEC, p. The freedom presented in the 
choice of directions of momenta results in a nonzero width 
of the otherwise (5-like structure of C2{6e) for ctq = 0. It 
then further broadens when one allows for some nonzero 
width (To- This means that it can be used as an addi- 
tional parameter when comparison with data would be 
attempted (provided its physical meaning will be made 
clear, see, for example, the discussion in Appendix [B] and 

El). 

If one wants to further continue with C2 [Sp^ ^ ^ ) , some 
additional input (new parameter(s)) is necessary, which 
has to be justified. We shall not discuss this point in 
detail, as it would bring us outside the main scope of 
this paper. Instead, we limit ourselves to showing in Fig. 
[7] the results of some more refined choices of directions 
of momenta. All results are for ctq = 0, introduction 
of nonzero ctq will change them accordingly in the man- 
ner presented in Fig. [51 They were obtained by choos- 
ing first values of transverse momenta, pt = \pt\, by 
selecting them from some exponential distribution con- 
strained only by the assumed mean value, (pt), which 
serves therefore as new parameter. This corresponds 
to selection of polar angles from the band centered on 

©Tjieari — arctan (^{pr) /pi"'"^''^ ), the corresponding axial 
angles were chosen uniformly from the [0, 27r] range. In 
this way, one gets components Px and Py and longitudinal 
component Pz = PL = ~ Pt- Two natural situa- 

tions are considered here: (i) - the maximally isotropic 
case (the pl of every consecutive particle, irrespectively 
of the EEC they belong to, is randomly assigned the sign 
(±)) and {ii) - the case in which bunching in energies is 
preserved also on the level of momenta (all particles in a 
given EEC have in the same hemisphere). All other 
choices should interpolate between these two. 




1 


* 5 isotropic 




5 " <p,j^>=3.o GeV 


a 


8 <p^>=o.3GeV 


\ 


P=o.9*e T=3.5 GeV 


1^ 

{(9 




A 




■ 








% 













1, 



* 5 isotropic 

— *— S <p^.>=3.o GeV 

6 <p.j,>=o.3 GeV 

P=o.9*e^'^, T=3.5 GeV 




[GeV] 



0.6 
0.0 



> 1 


9 5 isotropic 




—if— 5 ' <p.j,>=3.o GcV 


i 


5 <p.j,>=o.3 GeV 


4 


P^o.g'e"'^, T=3.5 GeV 



























[GeV] 



[GeV] 



FIG. 7: Comparison of two different ways of choosing mo- 
menta of particles occupying given EEC (|p| is always fixed) 
for (left panels) and 5p^ (right panels). Upper panels: 
the (+, — ) sign of p_L was chosen randomly for every particle 
without referring to EEC it belongs to. Lower panels: it was 
chosen randomly for EECs and kept the same for all particles 
in it. In all cases three choices of px is shown: unlimited 
[isotropic) and with limits imposed by two different values of 
(pt) when sampling pt from exponential distribution; in all 
cases (To = . 

Two characteristic features seen in Fig. [7]should be no- 
ticed: (z) one observes narrowing of C2{5p^) (i.e., trans- 
verse) distributions with tightening the allowed pT region 
(i.e., when proceeding from fully isotropic distributions 
to those restricted by assumed {pr) with diminishing val- 
ues of {pt}), this effect is essentially independent of the 
way the signs of Pz components are chosen, (ii) C2{Sp^) 
shows different behavior depending on which choice of 
signs is followed: it shows no dependence on (pt) for the 
choice (i) above whereas for the choice (ii) a difference 
shows up only for 6p^ > 0.2 GeV. The situation changes 
dramatically when one allows for smearing of energy in 
the EEC, i.e., for o-q > 0, see Fig. [8] This means intro- 
ducing a new parameter, to which our results are most 
sensitive but which is still not well understood (see, for 
example Appendix [B)) . 

In Fig. [9]results for C2{Se) and C2{Sp^) are compared 
with those for C2{Qinv), where Qi„D = - (Pi - ^2)^ with 
Pi_2 being 4-momenta of particles 1 and 2 (it is relativis- 
tically invariant variable , such that = implies for 
massive particles that pi = p2 and BEC has its maxi- 
mum there). Notice the peaked shape of C'2{Qinv) and 
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FIG. 8: The same as in Fig. [7] but for choice of pr restricted 
by (pt) =0.3 GeV and for different values of ao — 0.0, 0.1 
and 0.3. 
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FIG. 9: Comparison of C'2{5i = Se.p^) as presented in Fig. |6] 
for (To = 0.1 with C2{Si = Qinv) for the same parameters (left 
panel). Notice that maximum of C2{Qinv = 0) > 2 and to 
get it below this value one has to increase uo (see right panel 
with (To = 0.3). 

the fact that C2{Qinv = 0) > 2 here. As shown in [i^l 
both features are mostly due to the specific kinematics 
of the Qinv variable (because of which it collects in the 
first bin contributions from the whole range of momenta 
provide(i only that they are near enough to each other). 
The point is that in case of Qinv variable smearing a in 
energy E and smearing in momentum p are not indepen- 
dent and therefore they do not sum up (as is the case for 
independent variables) but rather they tend to cancel. In 
fact, their global effect is of the order of (J■{^/p^ + rn^—p) 



which tends to zero for large momenta and never exceeds 
a ■ m, where m is pion mass. In other words, smearing 
in Qinv is considerably smaller then that in energy and 
momentum. Therefore there is a natural tendency of in- 
creasing occupancy of the lowest bins in Si = Qinv in 
comparison to (5_e or Sp^ . 




5^ [GeV] 

FIG. 10: Comparison of C2(SE) for different maximally al- 
lowed sizes of EEC given by maximal number of particles 
n,nax they can have. Notice difference from the similar re- 
sults obtained for pionic lattice and presented in Fig. [2};. It 
is caused by the fact that here distance between particles is 
not limited by the fixed spacing of the lattice as before. 

Let us now stress the most specific feature of our al- 
gorithm: it always produces BEC of the highest possible 
order n (for which C2(Se = 0) = nl). This order is 
dictated by the the number of particles in the most pop- 
ulated EEC, which in turn depends strongly on location 
of a given EEC in phase space, see Table |TT1 It will signifi- 
cantly influence both the strength (as given by C2{S = 0)) 
and shape of the BEC effect, cf., Fig. [TOl In it we show 
what happens when the maximal number of particles in 
each EEC is artificially limited so as not to exceed some 
imposed maximal value equal to Umax- Notice that this 
requirement also affects the resulting number of EECs 
(the smaller population of EECs the more of them must 
be present). It is clear that this effect will be stron gest 
in events with very high multiplicities recorded (cf. |44l | 
for references to projects of the respective experiments). 
The sensitivity of our algorithm to rimax^s presented in 
Fig. [TO] (for C2{5p^) the changes are qualitatively the 
same for all choices of momenta presented here) makes 
it an ideal tool for numerical investigations of BEC also 
for particles satisfying statistics different than BE (for 
example, the so called parastatistics ^45]). 

We close this Section showing how sensitive are C2{Si) 
to different choices of EECs represented by different val- 
ues of parameters Vq and T and to different masses W of 
hadronizing source, cf. Fig. [TT] and Table [J (this is done 
again for (Tq = and using isotropic distributions of di- 
rections of momenta. Any changes in them, as discussed 
before, would then change these results accordingly in the 
way demonstrated in Figs. [HI [5] and [7]) . The most char- 
acteristic feature is the observed growth of C2{Si = 0) 
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FIG. 11: Upper panels: sensitivity of C2{Sp^) to different 
choices of EECs exemplifying by different sets of parameters 
used. Lower panel: the same but for different energies W 
and fixed size of EEC (as indicated, in all cases ao = and 
isotropic source was considered). The corresponding values 
of mean multiplicities, their dispersion as well as the mean 
number of EECs and their mean occupancies are listed in 
Table I] 

(i.e., the grow of the so called "parameter of chaoticity" 
X = C2{Si — 0) — I ) with diminishing number of EECs, 
{Nceii)- Actually, this result was behind the original in- 
troduction of the notion of EECs in describing the BEC 
effect done in [12] where the number of EECs is tightly 
connected with the parameter k in the NB multiparticle 
distributions used to fit data. Interesting feature of this 
approach, shared by our picture as well, is that, as shown 
in [12], it explains in a natural way the dependence of A 
on dN/dy and on the atomic number of projectiles A [46j . 
In Appendix [Cl we discuss changes in C2{Si) introduced 
when correcting for the inevitable energy-momentum and 
charge imbalances induced by our algorithm. 



C. Discussion 

Let us recapitulate the physical picture we are propos- 
ing. Its basic object is an EEC, a quantum state con- 
taining a number of identical secondaries of the same, 
or nearly the same energies. These secondaries are as- 
sumed to satisfy the BE statistics which is imposed by 
demanding that they follow geometric distribution. As 
a result the correlation functions C2{Sx) that follow are 



very sensitive to the characteristics of the EEC, for ex- 
ample the width of C2{Sx) is proportional to the allowed 
energy spread in EEC, a. It is best seen on the example 
of C2{Se), cf. Fig. [HI Interestingly enough, when all par- 
ticles in the EEC have the same energy, C2{Se), is diver- 
gent. On the other hand, in the same situation, correla- 
tion functions in momenta have already nonzero widths. 
This is because the choice ct = does not constraint di- 
rections of momenta. The simplest case of isotropically 
selected directions, represented by C2{Sp^), is shown in 
Fig. [51 The choice of directions provides therefore addi- 
tional freedom in modeling C2{Sp^ ^). It can be seen in 
Figs. [7] and [8] where examples of restricting the range of 
transverse momenta and the choice of longitudinal mo- 
menta are displayed. As for the physical meaning of 
CT > 0, it is tempting to identify it with the temporal 
characteristic of hadronization process, in fact with the 
life time of an average EEC. 

So far we presented results only for direct pions being 
produced. To include resonances one would first have to 
decide whether the BEC should affect them in the same 
way as particles or whether they affect only the pions 
resonances decay into. This point surely deserves further 
discussion, but it would bring us outside the scope of our 
paper. Therefore we present in Fig. [T2] results for C2 
calculated for p mesons (of charges (-1-,— ,0) and mass 
nip — 0.7 GeV, with zero widths assumed for simplicity) 
considered to be simple particles subjected to the same 
procedure of building EEC as that for pions. This is com- 
pared with the case where such p's subsequently decay 
into pions. Notice that whereas the BEC for p's is quite 
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FIG. 12: Comparison of BEC for mesons p (with rup — 0.7 
GeV) and for pions obtained from their decays. To maximize 
effect we use the EEC's for p's with fixed energies (i.e., ctq = 
0). Both C2{5e) and C2((Sp^) are presented. Notice that also 
for p correlation function C2{Se) is concentrated in first bin 
only. 

strong and not very much different from that for pions 
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only, it hardly survives the process of p's decay, espe- 
cially for the C2{Sp^) case. However, it must be stressed 
at this point that, so far, our p's were treated as spinless 
particles and that they were assumed to decay into pair s 
of pions in an isotropic way in their center of mass [47[ . 

In our algorithm only particles in EEC are subjected to 
BEC, there are no intercorrelations of the BE type be- 
tween particles from different EECs (see Appendix [X| . 
Such picture seems to be supported by recent data on 
BEC in e+e^ — > W^W^ multi-Vt^ boson production 
which show non existence of inter-iy BEC. This result 
suggest strongly that although spatially located practi- 
cally on top of each other, nevertheless bosons W act as 
independent sources of pions in this case [48| . 

We would like to stress here that all restrictions on 
energies and momenta mentioned above influence first of 
all the shape of a single EEC in momentum space rather 
than global characteristic of a given hadronizing source 
used. It is therefore to a large extent the EEC infor- 
mation on which are encoded in the correlation function 
C2 ■ On the other hand their number depends in the first 
instance on the characteristics of the hadronizing source 
encoded in the choice of f{E). Occupancy of EEC, how- 
ever, depends strongly on Vq and together with T from 
f{E) change substantially both the multiplicity distribu- 
tions (in our case from the Poisson-like to Polya-Aeppli 
one) and single particle distributions as well (see Figs. [J] 
andO. 
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M = Y/Sy M = Y/6y 

FIG. 13: Example of intermittency signals obtained using 
our algorithm with Vo = 0.9, T = 3.5. Left panel: not cor- 
rected for energy-momentum and charge imbalances results 
for (To = (black curves) and (tq = 0.3 (grey curves) are com- 
pared. Right panel: results for cro = only with and without 
corrections (most of the visible effect comes from correcting 
the energy momentum imbalance). Moments Fq{M) are de- 
fined as in 29]. 

It is worth to mention at this point that together with 
BEC our algorithm introduces also some intermittency 
signal, much in the way expected already in [2^, see Fig. 
[T3I It depends noticeably on the smearing of EEC in 
energy and is very sensitive to the energy momentum 
imbalance corrections (cf., Appendix [Cjl ^49l] . 



IV. SUMMARY AND CONCLUSIONS 

To summarize: we argue that proper numerical sym- 
metrization of the multiparticle state of identical par- 
ticles can be achieved in an economical way (in what 
concerns computational time) only by bunching them 
in phase space in such way that (identical) particles 
in each bunch (called here EEC - elementary emitting 
cell) have (almost) the same energies and follow a ge- 
ometrical (Bose-Einstein) distribution. Only particles 
in EEC experience BEC, those from different cells do 
not (see Appendix [A|. We regard this conjecture as 
emerging in a natural way from previous investigations 
[1, O O, [H, [131 . It is the main result of studies of 
a properly symmetrized multiparticle wave function of 
identical secondaries produced in the reaction [l^, [l^ . 
They unravelled that in such state the originally uni- 
formly distributed particles start to bunch [3, [l^, also 
changing the original poissonian multiplicity distribu- 
tion to a NB one The similarity with the clan 
model \2§\ leading to QCM proposed here was then im- 
mediate. The same conjecture was achieved indepen- 
dently following studies in which one works with num- 
ber of q uant a (p articles) without invoking wave functions 
[H [Ml m, [11 m m, 113. in the first practical appli- 
cation of the concept of bunching [l^ the whole (one 
dimensional) phase space was divided into such EECs 
(of equal size in rapidity space) and this simple decision 
resulted in profound consequences for what concerns the 
ability to describe different physical distributions. We 
develop this idea further: our EECs are formed dynam- 
ically, they can both overlap each other and be widely 
separated from each other and their number and multi- 
plicities (i.e., their sizes) of particles in them fluctuate 
from event to event. This work is then about how to 
form such EECs and what to do with them. 

One can ask whether cells (or proper bunching of iden- 
tical particles in phase-space, advocated here) are really 
needed to model quantum Bose statistics or, perhaps, 
is would be enough just to use the usual Bose-Einstein 
distributions to this aim. To answer this question let 
us compare two ways of producing bosonic particles: (i) 
generating them directly from Bose-Einstein distribution 
(like {N{Ei)) as given in eq. p^ . which is the most sim- 
ple way advocated on many occasions) and, (ii) generat- 
ing particles from a Boltzmann distribution and bunch- 
ing them in an appropriate way in phase-space according 
to QCM presented here. In the first case, one gets the 
correct single particle distribution, whereas in the sec- 
ond case one also accounts for multiparticle correlations 
of quantum statistical origin. The corresponding results 
are presented in Fig[T31 They show that both approaches 
lead to very different results, which can be understood by 
realizing that case (z) corresponds to a particular realiza- 
tion of case (ii), namely with only one cell containing the 
same average number of particles. In fact case (i) shows 
only some trivial correlations which can be eliminated by 
a proper choice of the reference distribution. Therefore, 
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FIG. 14: Comparison of C2 modelled by using MC event gen- 
erators with EEC (circles) and without EEC's but with select- 
ing particles directly from the corresponding Bose-Einstein 
distribution {N{Ei)) as given in eq. (I14|> [soj . In both cases 
the Boltzmann distribution for classical particles was used as 
reference event. This fact may be crucial to get apparent in- 
crease in case of using directly Bose-Einstein distributions (as 
in [50|). The use of mixed event instead, which will not af- 
fect substantially the result obtained using EECs, will most 
probably kill the other effect. 



we argue that, according to our understanding, the bet- 
ter approach (after correcting additionally for charge and 
energy-momentum conservations) is that 



BEC = cells + geometrical distribution 



(20) 



From its construction it is evident that our algorithm is 
best suited to study events with large multiplicities (as 
planned in some experiments [11] )■ The statement ([^01) 
also summarizes the only possible way to incorporate our 
algorithm in MCEG codes: according to our finding it 
should be done by enforcing particles to be produced in 
bunches with characteristics of our EECs. In fact, closer 
inspections on all previous efforts to imitate BEC men- 
tioned in Section Hi Al shows that they all implicitly were 
aiming in similar direction [13, El, 0.9, 20, 21]. The dif- 
ference was that they were trying to do it a posteriori 
rather than a priori and that their weighting procedure 
was aimed more at the spatio-temporal (unknown) char- 
acteristics of the hadronizing source than on the true 
physical principles of BEC. 

In this work we have stopped short of a general com- 
parison with data. The reason is that we were using here 
a most simple statistical model of hadronization in order 
to illustrate our algorithm. According to it a hadronizing 
source is assumed to be a single object with some fixed 
mass W and fixed initial charge, which is the situation 
encountered only in e~^e~ annihilation reactions. In all 
other multiparticle production processes one either en- 
counters W varying from event to event (foUowing^ome 
distribution, for example, inelasticity distribution 51]) or 
there is a number of hadronizing sources with different 



W in each event (which is the most probable situation in 
heavy ion collisions). Moreover, the statistical scheme of 
hadronization employed here means that our hadroniz- 
ing source does not experience any internal flows and is 
not subjected to any external force, which would result 
in some energy-position correlations or specific effects of 
partial coherence (cf., [l^) and thus additionally influ- 
ence the C2 correlation function. The problem of the 
possible net charges of such sub-sources was never dis- 
cussed. All this asks for very specialized study, which 
goes outside the scope of this paper. What we can only 
say at this moment is that, when undertaken, such a 
study would mean the necessity of developing our formal- 
ism further in what concerns details of modeling EECs 
mentioned before. The most probable approach would 
be to additionally assume that each EEC itself should be 
described by a properly symmetrized ripart-particle wave 
function (where Upart is the number of particles in this 
EEC). Only then would one be able to introduce into 
the model a characteristic momentum-positions correla- 
tions caused by quantum statistics (much in the spirit 
leading to eq. ([T])). The price for following such a proce- 
dure is the necessity to also introduce to our description 
the space-time characteristics of the hadronizing source, 
which we were not dealing with so far. This would re- 
sult in the same permanent structure as given by eq. 
but with much smaller sizes and with explicit dependence 
on spatio-temporal variables (as in eq. ([8])) to be used 
later when selecting momenta (actually, they would be 
effectively integrated out in this procedure). However, 
the noticeable feature is that in our case this procedure 
would not demand spatio-temporal factorization prop- 
erty of the hadronizing source assumed in ([1]). When 
replacing plane waves used there by Coulomb distorted 
wave functions |52i] (which is common practice nowadays) 
one could then attempt, in principle, to account for the 
influence of Coulomb interactions so far neglected here 
(albeit only on the level of 2-body interactions). This is, 
of course, not the only possible generalization and there- 
fore we leave the problem of confrontation with real data 
for further studies. Finally, let us notice that any seri- 
ous comparison with data would have to be done includ- 
ing corrections for energy-momentum and total charge 
imbalances induced by our algorithm, which (as seen in 
Appendix [C| can be quite substantial and not unique. 

Let us close by noticing that, although our algorithm 
was originally intended to model quantum phenomena 
of BEC only, it is in fact more general. The reason is 
that the characteristic structure of the C2 correlation 
function associated with specific bunching of identical 
particles turns out to be quite universal phenomenon 
also observed in many other, purely classical systems, 
provided only that they exhibit strong and correlated 
fiuctuations [53] • This means then that our algorithm, 
albeit with different motivation, could also be applied 
there. At this point one should also notice an attempt at 
numerical modeling of another quantum phenomenon, 
namely Bose-Einstein condensation presented in [s^ ] and 
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the possible connection of the above with the physics of 
networks |55|]. Finally, we also claim that, because of its 
sensitivity to maximal occupancy of EECs, out method 
could easily be modified to be able to study BEC effects 
for para-bosons [isj . 



as degree of coherence. On the other hand we can write 
following Sec. Ill CI the true correlation function C2 as 



n{n-l) 

O2 — — i 



Var{n) 



1 

n' 



(A5) 
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APPENDIX A: SOME REMARKS ON EEC 

We would like to comment in more detail on our propo- 
sition that identical bosons should come in EEC and 
experience effects of BEC there, whereas no such ef- 
fects should be seen when bosons from different cells are 
considered. In the language used in [56] the probabil- 
ity of registration of a coincidence of n bosons in states 
ji , . . . , jn is given by a normalized correlation tensor of 
rank 2n, 



(nn) 

9jl--.jnjn+l...j2n 



Q{nn) 

n-3^3^^.-n^ where 



n 



In 



G 



(11) 



(Al) 



G 



(nrn) 



Tr 



is given in terms of the density matrix operator p whereas 
and a are, respectively, creation and annihilation op- 
erators. If 



1 



n< N 



(A2) 



we have coherence of the order 2N . Experimentally this 
means that the probability of registering n bosons in co- 
incidence is equal to the product of probabilities to reg- 
ister individual ones. Because of commutation relations 



for bosons, [ak, ai] 

has that a\a\akai 
n^n/ and then from eq 
different states 

(22) 



and 



ak,al 



Ski, one 



ajalakai = a\a\aiak = a\a\aiak = 
ll) that for two bosons from 



rikni 
nk ■ W 



rik ■ ni 
nk ■ ni 



= 1. 



(A3) 



where Var{n) = Var{nk) = n(l -I- n) for bosons from 
the same state, which in our case means from the same 
EEC. This immediately leads to C2 = 2 in this case (for 
Boltzmann statistics Var{n) — n and one gets 6*2 = 1 
instead). For k EECs we have multiplicity distribution in 
NB form with variance Var{n) = kVar(nk) — n -\- r? jk 
and mean multiplicity n = kfik- This leads immediately 
to C2 = 1 + 1/fc, i.e., we have C2 = 2 for /c = 1 and 
C2 — > 1 for k — » 00. 



APPENDIX B: STRUCTURE OF C2(Q) AND 
FINITENESS OF HADRONIZING SOURCE 

As demonstrated in [lol | using a quantum field the- 
ory approach to BEC, the structure of the correlation 
function G2{Q) is connected with the finiteness of the 
hadronizing source. In the four-dimensional case, the 
wave function formalism with exp (—ip^ ■ x'^) implies 
the infinite (4-dimcnsional) volume of the hadronizing 
source, V 00. This can then be connected to the fact 
that commutation relations for the respective operators 
contain delta functions, 



(Bl) 



which are nonzero only for identical values of the four- 
momenta (i.e., also energies), . They in turn lead to 
correlation functions C2{Q) in the form of 



C2(Q) = 1 + 5(0 -i?). 



(B2) 



i.e., to C2{Q) > 1 but only at one point, Q = (cf., 
Table HlB . However, as shown in [IQ] . assuming commu- 
tation relations in the form of some sharply piked (but 
not infinite) functions A, i.e., replacing delta functions 
by functions with supports larger than limited to a one 
point only. 



This means that two bosons from different states (in our 
case: from different EECs) exhibit second order coherene. 
On the other hand, in the situation in which only one 
state k is occupied, i.e., when alaluka 
2 o„ iIaTI) results in 



alak 



alttkalak 



nk, eq. 



(22) _ nf. - Uk _ 2nk _ 
9k,k,k,k, - — — - _._ " 



(A4) 



(here we use the fact that in geometrical distribution nf, — 
TTk^ = nfe(l + rik) ). Notice that it is greater by unity 
than (jA3|) but, because g is not limited, it cannot serve 



[c{p^),c{pI)]^AHp^-pI), 



(B3) 



results immediately in the correlation function endowed 
with final width: 



C2(0) - i + fl(0-i?) 



(B4) 



where the form of g((5 • i? is a straightforward reflection 
of the assumed form of A''(p^ — p'^)). Such a proce- 
dure corresponds to introducing a finite dimension of the 
source, V — Vq (and the use of wave packets formalism, 
exp [—ipfj. ■ — p^ / (2o'p)] , instead of plane waves). 



TABLE III: Schematic presentation of how a different form of 
commutation relations in a quantum field theoretical descrip- 
tion of BEC results in a different form of C2 function jlO| |. 
The actual form of C2 (i-e., of the function • R)) depends 
on the form of function A(. . . ) used to moderate the original 
(5(. . . ) function. 



Volume 
(4dim) 
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Correlation 
function 
C2(Q) - 1 


V ^ 00 
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e ^"P 
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APPENDIX C: CORRECTING FOR 
ENERGY-MOMENTUM AND CHARGE 
IMBALANCES 

Any random selection procedure, even when follow- 
ing formulas which assume exact energy-momentum and 
charge conservations, frequently induces some energy- 
momentum and charge imbalances and this is also true in 
the case of our algorithm |52| . The obvious remedy would 
be to only accept events preserving the initial values of 
energy-momentum and charge. This, however, would re- 
sult in an unacceptable long computational time. To the 
best of our knowledge, only in the algorithm presented 
in Section III B 31 [16] is energy-momentum assured us- 
ing this method (with < 20% accuracy) and charges are 
conserved by accepting only events reproducing the ini- 
tially assumed charge. Results of all other algorithms 
presented here were prepared in the same way as what 
we have presented, i.e., not corrected for any, therefore 
they can be compared with each other. We would now 
like to discuss changes in correlation functions induced 
by correcting for energy-momentum and charge imbal- 
ances introduced in the selection process. At first one 
must stress that there is no unique procedure to per- 
form such corrections [SS^. In what follows we shall use 
a very simple (but fast) procedure: (z) - shifting (by the 
same amounts, A.px, Apy and Ap^) components of mo- 
menta, Px,y,z and, after balancing momenta, appropri- 
ately rescaling energies; (m) by converting the necessary 
number of (+) particles to (— ) ones (or vice versa, de- 
pending on the actual charge balance in the event) and, 
in case of odd number of particles with surplus charge, 
attribute the surplus charge to some randomly chosen (0) 
charged particle. 

The results are shown in Fig. \TE\ (our hadronizing 
sources have zero initial charge). The most sensitive to 
the correcting procedure is C2(fe), correlations of mo- 
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FIG. 15: Effects of correcting for the imbalances introduced 
by our selection procedure in energy-momentum (upper pan- 
els and lower-right panel) and total charge (lower panels). 
Results without (w/o) corrections are the same at all pan- 
els. Two selection of the (-f / — ) sign of the pl component are 
used: it is chosen randomly for each particle irrespective to 
EEC it belongs to (all panels except the upper-left one); all 
particles in a given EEC possess pl of the same sign (upper- 
left panel). In all panels first bin of gray full curve contains 
uncorrected C2[5e)- In all cases cro = 0. 




FIG. 16: Typical distributions of total charges obtained from 
our selection procedure performed for W = 91.2 GeV, T = 3.5 
GeV and cro = 0. for different sizes of EEC as dictated by 
parameter Vo — 0.9 and Vo = 0.3 and for maximal size of EEC 
artificially limited to Umax = 4 (for Vo — 0.9). The best fits 
to most narrow and most wide distributions are shown. Total 
initial charge of hadronizing mass W was assumed Q = 0. 



menta feel corrections only when all particles in a given 
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EEC are located in the same hemisphere. Otherwise 
there is essentially no difference, The reason for such be- 
havior is that in the other case, the relative differences of 
momenta considered here are not changed by the shift- 
ing procedure. In correcting for AQ 7^ the EECs with 
single particle only were chosen first, afterwards EECs 
with Npart > 1 were chosen randomly until the correct 
balance of charge was achieved. In this case, as one can 
see in Fig. [15] (lower panels), the effect is quite dramatic 
and essentially independent of the way the momenta were 
chosen or on their energy-momentum balance. The best 
measure of this is provided by the widening of the origi- 
nally S-\ike 02(6 e)- Notice also the clearly visible upward 
bending of the tail of C2 distributions. It is worth notice 
at this point that similar shapes are observed in C2 ob- 
tained from e+e~ annihilation experiments, in which, as 
in the case considered here, the original energy and total 
charges are well known and fixed. 

Effects shown here are so dramatic because, in order 
to be as near as possible to the proper BE distributions 
in the average EEC in situation of only limited energy 
W available for hadronization, we had maximized sizes 
of EECs by allocating to them many particles, this was 
technically achieved by using large value of parameter 
'Po, 'Po — 0.9. It resulted in very broad total charge dis- 
tribution (centered on the assumed value Qtot — 0), 



see Fig. [161 The reason for such large fluctuations is 
as follows. In our algorithm each EEC contains parti- 
cles of the same charge (-I-, — , 0) selected randomly (with 
equal probabilities). With only one particle per EEC, 
Npart — 1, this would result in quite narrow AQeec- 
However, because in general Npart > 1 and fluctuates, 
the charge imbalance AQeec broadens considerably to 
what is observed in Fig. 1161 It leads then to very large 
fluctuations and events with large charge imbalance are 
quite frequent. The effect is therefore, as seen in Fig. \W[ 
very sensitive to the size of EEC allowed (see Fig. [Hand 
Table [H, i.e., to the parameter Vq responsible for Npart 
and to any attempts to limit it as, for example, shown in 
Fig. [10] [HI . 

To summarize: this problem arises because our pro- 
cedure destroys part of the originally formed EECs and 
forms some new ones what results in the sensitivity ob- 
served. On the other hand we do not see at the moment 
any economical way to account for extremely complicated 
correlations arising when attempting to keep Qtot = 
all the time. The only apparent cure, to keep only events 
with exactly right charge balance, would place our algo- 
rithm in the same category (in what concerns the use of 
computational time) as those presented in Sects. IIIB II 
and lIIB2] 8l.[l5l| (which, by the way, were not attempting 
to impose any conservation laws at all). 
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